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Abstract 

A  preliminary  algorithm  based  on  the  two-fluid  plasma  model  is  developed  to  in¬ 
vestigate  the  possibility  of  simulating  plasmas  with  a  more  physically  accurate  model 
than  the  MHD  (magnetohydrodynamic)  model.  The  algorithm  is  based  on  a  Roe-type 
approximate  Riemann  solver.  Beginning  with  the  two-fluid  plasma  model,  the  govern¬ 
ing  equations  are  normalized  and  formulated  in  conservation  form.  The  eigenvalues  and 
eigenvectors  of  the  system  flux  Jacobians  are  determined  and  properly  normalized  to 
prevent  catastrophic  cancellation.  An  approximate  Riemann  solver  is  developed  based 
on  the  derived  conserved  fluxes.  The  electromagnetic  fields  are  solved  and  coupled  to 
the  two-fluid  approximate  Riemann  solver.  The  field  solver  is  an  electrostatic  potential 
solver  appropriate  for  electric  fields  and  solves  Poisson’s  equation.  The  algorithm  is 
benchmarked  against  analytical  results  and  published  simulations.  In  a  parallel  effort, 
numerical  methods  for  the  calculation  of  the  flux  Jacobian  were  be  investigated.  The 
primary  methods  are  the  first-order  limit  formulation  and  the  second-order  complex 
number  formulation. 
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1  Executive  Summary 

The  primary  objective  of  this  project  is  perform  the  preliminary  work  to  develop  a 
plasma  simulation  algorithm  based  on  the  two-fluid  plasma  model.  A  two-fluid  plasma 
simulation  algorithm  would  be  capable  of  time-dependent,  three-dimensional  simula¬ 
tions  of  plasma  phenomena  in  realistic  geometries.  The  algorithm  would  be  a  valuable 
tool  for  the  design  and  testing  of  plasma  related  technologies  that  are  important  to 
the  Air  Force  and  industry,  such  as  advanced  plasma  thrusters  for  space  propulsion, 
hypersonic  drag  reduction,  portable  pulsed  power,  high  power  microwave  devices,  nu¬ 
clear  weapons  effects  simulations,  radiation  production  for  counter  proliferation,  and 
fusion  for  power  generation.  Implementing  the  algorithm  on  parallel  supercomputers 
would  allow  the  detailed  modeling  of  realistic  plasmas  in  complex  three-dimensional 
geometries. 

Current  plasma  simulation  codes  capable  of  complex  geometries  are  based  on  the 
MHD  (magnetohydrodynamic)  model.  The  derivation  of  the  MHD  model  involves  sev¬ 
eral  assumptions  that  severely  limit  its  applicability.  The  two-fluid  model  only  assumes 
local  thermodynamic  equilibrium.  A  two-fluid  algorithm  is  more  physically  accurate 
and  capable  than  MHD  models.  The  two-fluid  model  is  normalized  and  formulated  in 
conservation  form.  The  algorithm  uses  an  approximate  Riemann  solver  that  we  devel¬ 
oped  for  the  two-fluid  plasma  model.  The  algorithm  is  benchmarked  against  known 
analytical  results  and  published  simulations. 

Several  professional  collaborations  exist  between  the  Department  of  Aeronautics  and 
Astronautics  at  the  University  of  Washington  and  the  Air  Force  Research  Laboratory 
both  at  Kirtland  AFB  and  Edwards  AFB,  Los  Alamos  National  Laboratory,  Lawrence 
Livermore  National  Laboratory,  the  University  of  Michigan,  Stanford  University,  and 
other  departments  at  the  University  of  Washington.  The  work  from  this  project  has 
been  presented  at  international  conferences. 


2  Project  Description 

Plasmas  are  essential  to  many  technologies  that  are  important  to  the  Air  Force,  some 
of  which  have  dual-use  potential.  These  applications  include  portable  pulsed  power  sys¬ 
tems,  high  power  microwave  devices,  drag  reduction  for  hypersonic  vehicles,  advanced 
plasma  thrusters  for  space  propulsion,  nuclear  weapons  effects  simulations,  radiation 
production  for  counter  proliferation,  and  fusion  for  power  generation.  In  general,  plas¬ 
mas  fall  into  a  density  regime  where  they  exhibit  both  collective  (fluid)  behavior  and 
individual  (particle)  behavior.  The  intermediate  regime  complicates  the  modeling  of 
plasmas. 


2.1  State  of  the  Art  —  Kinetic,  PIC,  and  MHD  Models 

Plasmas  may  be  most  accurately  modeled  using  kinetic  theory.  The  plasma  is  described 
by  distribution  functions  in  physical  space,  velocity  space,  and  time,  /(x,  v,t).  The 
evolution  of  the  plasma  is  then  modeled  by  the  Boltzmann  equation. 
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for  each  plasma  species  a  =  ions,  electrons.  The  Boltzmann  equation  coupled  with 
Maxwell’s  equations  for  electromagnetic  fields  completely  describe  the  plasma  dynamics. 
[1-3] 

However,  the  Boltzmann  equation  is  seven  dimensional.  As  a  consequence  of  the 
large  dimensionality  plasma  simulations  using  the  Boltzmann  equation  are  only  used  in 
very  limited  applications  with  narrow  distributions,  small  spatial  extent,  and  short  time 
durations.  [4, 5]  The  seven  dimensional  space  is  further  exacerbated  by  the  high  velocity 
space  that  is  unused  except  for  tail  of  the  distribution  or  energetic  beams.  Boundary 
conditions  are  difficult  to  implement  in  kinetic  simulations. 

Particle  in  cell  (PIC)  plasma  model  apply  the  Boltzmann  equation  to  representative 
superparticles  which  are  far  fewer  (10^)  than  the  number  of  particles  in  the  actual 
plasma  (10^°).  [6]  PIC  simulations  have  similar  limitations  as  simulations  using  kinetic 
theory  due  to  statistical  errors  caused  by  the  fewer  super  par  tides.  Boundary  conditions 
are  also  difficult  to  implement  in  PIC  simulations. 

The  other  end  of  the  spectrum  in  plasma  model  involves  taking  moments  of  the 
Boltzmann  equation  and  averaging  over  velocity  space  for  each  species  which  implicitly 
assumes  local  thermodynamic  equilibrium.  The  resulting  equations  comprise  the  two- 
fluid  plasma  model.  The  two-fluid  equations  are  then  combined  to  form  the  MHD 
(magnetohydro dynamic)  model.  [7]  However,  in  the  process  several  approximations  are 
made  which  limit  the  applicability  of  the  MHD  model  to  low  frequency  and  ignores  the 
electron  mass  and  finite  Larmor  radius  effects. 

The  MHD  model  treats  the  plasma  like  a  conducting  fluid  and  assigning  macro¬ 
scopic  parameters  to  describe  its  particle-like  interactions.  Plasma  simulation  codes 
based  on  the  MHD  model  have  been  very  successful  in  modeling  plasma  dynamics  and 
other  phenomena.  Codes  such  as  MACH2  are  based  on  arbitrary  Lagrangian/Eulerian 
formulations.  [8]  ALE  codes  are  well  suited  for  simulating  plasma  phenomena  involving 
moving  interfaces.  [9]  However,  ALE  codes  cannot  be  formulated  as  conservation  laws 
and  lack  many  of  the  inherent  conservative  properties. 

We  have  successfully  implemented  the  MHD  model  in  conservative  form  to  simulate 
realistic  3-D  geometries.  [10-12] 


dQ 

dt 
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where  Q  is  the  vector  of  conserved  variables,  Th  is  the  tensor  of  hyperbolic  fluxes,  and 
Q Non-ideal  Contains  the  non-ideal  (mostly  parabolic)  terms.  (For  algorithm  details  we 
refer  the  reader  to  our  published  articles,  Ph.  D.  dissertations,  and  previous  AFOSR 
technical  reports.  [10-16])  We  formulated  the  algorithm  implicitly  and  demonstrated 
the  superior  accuracy  and  convergence  over  an  explicit  mode.  Figure  1  shows  the 
convergence  history  using  the  implicit  algorithm.  The  algorithm  uses  a  finite  volume 
implementation  with  a  block-structured  domain  that  accommodates  arbitrary  shapes. 

The  lack  of  a  divergence- free  magnetic  field  is  one  of  the  most  severe  difficulties  in 
algorithms  which  use  the  MHD  model.  The  divergence  of  the  magnetic  field  must  be 
identically  zero  according  to  Maxwell’s  equation,  V  •  B  =  0.  A  common  remedy  is  to 
perform  a  Hodge  projection  to  “clean”  the  divergence  component  of  the  field.  We  have 
added  a  divergence  cleaner  to  our  MHD  algorithm.  However,  this  adds  a  nonphysical 
step  which  still  does  not  completely  eliminate  the  divergence. 

An  even  more  severe  limitation  of  the  MHD  model  is  the  treatment  the  Hall  effect 
and  diamagnetic  terms.  These  terms  represent  the  separate  motions  of  the  ions  and 
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Figure  1:  Convergence  history  using  the  SGS  method  to  invert  the  implicit  operator,  n  is 
the  number  of  physical  time  iterations,  m  is  the  number  of  pseudo  time  subiterations,  and 
sgs  is  the  number  of  iterations  of  the  SGS  method. 
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Figure  2:  Dispersion  relation  for  the  fast  and  slow  modes  with  uJcita  =  3.5.  The  dashed 
lines  represent  the  ideal  MHD  modes. 


electrons  that  still  preserve  the  plasma  approximation  of  quasineutrality.  The  Hall  effect 
and  diamagnetic  terms  also  account  for  ion  current  and  the  finite  Larmor  radius  of  the 
plasma  constituents.  These  effects  are  important  in  many  applications  such  as  electric 
space  propulsion  thrusters:  Hall  thrusters,  magnetoplasmadynamic  (MPD)  thrusters, 
Lorentz  force  thrusters.  The  Hall  term  is  also  believed  to  be  important  to  electrode 
effects  such  as  anode  and  cathode  fall  which  greatly  affect  many  directly  coupled  plasma 
applications. 

In  general  the  Hall  terms  are  difficult  to  stabilize  because  they  lead  to  the  whistler 
wave  branch  of  the  dispersion  relation.  The  phase  and  group  velocities  of  the  whistler 
wave  increase  with  frequency.  See  Figure  2  for  the  dispersion  with  a  modest  value  of 
the  Hall  parameter  (cx^ciTa)”^. 

We  have  treated  the  Hall  term  using  a  semi-implicit  technique.  [17, 18]  After  the 
hyperbolic  terms  of  eqn(2)  are  advanced,  the  Hall  terras  are  treated  independently. 
The  conserved  variables  are  then  corrected.  The  procedure  can  be  computationally 
intensive.  The  operator  stencil  uses  5  points  in  the  sweep  direction  and  3  points  in 
the  orthogonal  direction.  The  complete  operator  stencil  is  45  points.  We  have  found 
that  the  semi-implicit  method  works  adequately  for  small  Hall  parameters,  but  becomes 
unstable  or  slow  to  converge  for  large  Hall  parameters  often  seen  in  applications. 

We  have  added  multi- temperature  capabilities  to  our  MHD  code.  [16]  The  algorithm 
is  based  on  an  approximate  Riemann  solver  to  advance  the  conserved  variables.  The 
various  energy  loss  and  transfer  mechanisms  are  then  used  to  account  for  the  change 
in  temperature  in  each  of  the  constituent  species  of  the  plasma  and  neutral  gas.  The 
multiple  temperatures  are  then  used  to  modify  the  ionization  and  recombination  rate 
parameters.  However,  the  (primary)  hyperbolic  algorithm  is  still  a  single  fluid,  single 
temperature. 
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2.2  The  Next  Step:  Two-Fluid  Plasma  Model 

The  logical  next  step  to  extend  the  realistic  geometry  capabilities  of  the  MHD  model  to 
the  more  physically  accurate  two-fluid  model.  The  complexity  of  the  two- fluid  model  is 
greater  the  MHD  model  but  significantly  less  than  the  kinetic  model. 

The  two-fluid  model  is  derived  by  taking  moments  of  the  Boltzmann  equation  for 
each  species.  The  process  of  taking  moments  eliminates  the  velocity  space  and  yields  a 
representative  fluid  velocity  for  each  species.  The  model  has  the  same  dimensionality 
as  the  MHD  model  except  there  are  two  fluids.  (The  details  of  the  model  are  described 
in  Section  2.4.1.) 

The  physical  accuracy  of  the  two-fluid  model  is  much  better  than  the  MHD  model. 
The  only  approximation  made  is  local  thermodynamic  equilibrium  of  each  fluid  but  not 
with  each  other  fluid. 

The  two-fluid  plasma  equations  are  normalized  and  formulated  in  conservative  form. 
A  high  resolution  approximate  Riemann  solver  is  derived  that  accurately  calculates 
the  transport  of  the  conserved  quantities.  The  algorithm  is  implemented  using  finite 
volumes  to  allow  for  extending  the  algorithm  to  three-dimensions  and  to  model  realistic 
general  geometries  with  body  fitted  computational  grids. 

2.3  Research  Objectives 

The  objectives  of  the  proposed  project  are  to: 

•  Derive  an  approximate  Riemann  solver  for  the  two-fluid  plasma  model. 

•  Investigate  possible  methods  for  numerical  evaluation  of  the  flux  Jacobian. 

•  Determine  an  appropriate  flux  limiter  for  the  approximate  Riemann  solver  to 
achieve  a  higher  order  accuracy  algorithm. 

•  Formulate  an  electromagnetics  algorithm  and  couple  the  algorithm  to  the  approx¬ 
imate  Riemann  solver. 

•  Implement  the  algorithm  with  finite  volume  block-structured  grids. 

•  Validate  the  code  with  analytical  results. 

2.4  Technical  Description 

This  section  describes  the  details  of  the  two-fluid  plasma  model  beginning  from  the  two- 
fluid  equations  derived  from  moments  of  the  Boltzmann  equation.  Further  details  about 
the  derivation  of  the  model  can  be  found  in  most  advanced  plasma  physics  texts.  See  for 
example  Ref,  [19],  The  proposed  conservative  formulation  and  algorithm  are  presented. 
The  formulation  has  many  properties  which  simplify  the  nature  of  the  algorithm.  These 
properties  are  discussed. 

2.4.1  Two-Fluid  Plasma  Model 

The  two-fluid  plasma  model  is  derived  by  taking  moments  of  the  Boltzmann  equation 
[eqn(l)]  for  each  species.  The  process  of  taking  moments  averages  over  velocity  space 
to  yield  a  characteristic  fluid  velocity  for  the  distribution  of  each  species.  The  averaging 
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process  implicitly  assumes  local  thermodynamic  equilibrium  of  each  fluid.  The  moments 
are  taken  as: 


^  +  i^(E  +  v„xB). 
ax  rUa 


dv  = 


dU 

dt 


colliaions 


where  X  =  l,maVa,ma  (va  -  iy)af  for  each  plasma  species  a  =  ions,  electrons.  The 
fundamental  equations  for  the  two-fluid  model  are  generated  by  taking  moments  of  the 
Boltzmann  equation,  and  the  fundamental  variables  are  generated  by  taking  moments 
of  the  distribution  function. 

Maxwell’s  equations  describe  the  evolution  of  the  electromagnetic  fields.  The  fields 
couple  to  the  ion  and  electron  fluid  variables  and  complete  the  model. 

V  -  E  =  e(ni  -ne)/eo 
VB  =  0 


9B 

dt 
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=  V  xB//Xo-j 
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The  evolution  of  the  particle  density  of  the  ions  and  electrons  is  expressed  by  continu¬ 
ity  equations  with  the  source  term  equal  to  the  difference  of  electron  and  ion  ionizations 
and  recombination.  This  is  also  the  zeroth  moment  of  the  Boltzmann  equation. 


(5) 


=  nnni{(Tv)  i^ion  +  e-ion  ~  TlinQ{av)  recomb 


-^+V.(».v.).S.,.  (6) 

where  Tii,  Ue,  Tin  3-^6  the  ion,  electron,  neutral  number  densities,  v^,  Vg,  are  the  ion 
and  electron  fluid  velocities,  and  {p'v)  are  the  interaction  rate  parameters  for  ion  and 
electron  impact  ionization  and  recombination.  Partial  current  densities  can  be  defined 
as 

=  en^Vi,  and  je  =  engVe-  (7) 

where  e  is  the  electronic  charge.  The  derivation  presented  here  assumes  singly  ionized 
species,  q  =  e.  (For  multiple  charged  ions,  Ug  Zn^  where  Z  is  the  charge  state.) 
Using  the  partial  current  densities,  the  particle  continuity  equations  are  then  written 
as 


dfl-i  1  ^ 

(8) 

dt  e 

(9) 
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The  first  moment  of  the  Boltzmann  equation  yields  momentum  equations  for  each 
species, 

■ 

niTUi  -^  +  (vi-V)vi  = 

-  Vpi  +  TliC  (E  +  Vi  X  B)  +  Rei  -  niTtliViVin  (10) 

'5Ve  / 

neTUe  +  = 

-  Vpe  -  nee  (E  +  Ve  X  B)  -  Rei  “  neTTleVel/en  (H) 

where  E  and  B  are  the  electric  and  magnetic  fields,  pi  and  pe  are  the  ion  and  electron 
partial  pressures,  Rei  is  the  electron  to  ion  momentum  transfer  vector,  and  i/^n  and  i/en 
are  the  ion-neutral  and  electron-neutral  collision  frequencies. 

Substituting  the  definitions  of  partial  current  densities  and  using  eqns(8)  and  (9), 
the  momentum  equations  are  rewritten  in  conservative  form  for  the  partial  current 
densities. 


GTiif>  nig  / 


^E-— jeXB  +  — Rei+jeUen  +  -^  (13) 

\  Ti  ^  I 


The  second  moment  of  the  Boltzmann  equation  yields  energy  equations  for  each 
species. 

|ni  ^  +  (v»-V)T<  =-PiV-Vi-V -^i^-Pexu-Pc^  (14) 


—  Tie  ^  0^  +  (Ve  •  V)rej  —  — PeV  •  Vg  —  V  ■  qe  +  Pextc  Prad  (15) 

where  Ti,  Tg  are  the  ion  and  electron  temperatures,  qi,  qc  are  the  ion  and  electron  heat 
fluxes,  Pexti,  Pext^  are  external  input  powers,  Pcx  is  the  charge  exchange  power,  and 
Prad  is  the  radiated  power.  The  temperatures  are  related  to  the  partial  pressures  by 

Pa  =  n^Ta  (16) 

for  a  =  {i,  e}.  The  energy  equations  can  be  combined  with  the  corresponding  momen¬ 
tum  equations  to  yield  energy  equations  expressed  in  conservative  form  for  the  total 
energy. 

+  V  .  f(e.  +p,)  =  j,  ■  fE+  (17) 


■V)tJ  =- 
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dee 

dt 


V- 


(Ce  +  Pe)  — 

CTTg 


=  j..(E+^) 
\  euej 


f)  ^en  ~r 


where  the  total  energy  is  defined  by 


and 


3  „  1  2  3  1  2 

Si  =  -UiTi  +  -TiiTniVt  =  -Pi  +  -niTTiiVi 


3  1  2 

fe  =  -^Pe  +  -neTtleV^. 


(18) 


(19) 


(20) 


2.4.2  Conservative  Algorithm 

The  two-fluid  plasma  model  is  now  expressed  in  conservative  form. 

^  +  V  •  f  h  =  Qforce  (21) 

ot 

where  Q  is  the  vector  of  conserved  variables,  Th  is  the  tensor  of  hyperbolic  fluxes,  and 
Qforce  is  the  forcing  function.  The  vector  of  conserved  variables  is 
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(22) 


The  hyperbolic  flux  tensor  in  the  x  direction  is 
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(23) 


The  characteristic  velocities  will  then  be  calculated  to  construct  the  approximate  Rie- 
mann  fluxes. 
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In  1-D  with  an  assumed  ratio  of  specific  heats  7  =  5/3,  the  flux  Jacobian  for  the 
two-fluid  equations  is 
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The  eigenvalues  of  the  flux  Jacobian  give  the  characteristic  velocities. 

I  en,  ^  eni  y  3  ’  eng  ’  erig  Y  3  rrie  j 


(25) 


Acoustic  speeds  are  defined  as 


3  mi  3  miTii 


2  ^  ^ 

3  mg  3  meTie 


5  Pe 


(26) 


The  presence  of  a  complete  set  of  eigenvalues  proves  that  the  system  is  hyperbolic.  An 
approximate  Riemann  solver  is  appropriate.  The  eigenvalues  for  the  two-fluid  plasma 
model  represent  the  combination  of  the  drift  speeds  and  thermal  speeds  for  the  electrons 
and  ions.  The  speeds  for  the  electrons  are  expected  to  be  quite  large,  but  much  less  than 
the  speed  of  light.  The  large  span  of  the  eigenvalues  makes  the  system  stiff  and  more 
difficult  to  solve,  particularly  with  an  implicit  method.  The  stiffness  of  the  equations, 
particularly  the  source  terms,  is  solved  by  treating  the  source  terms  implicitly.  The 
details  are  described  in  Sec.  3.1.3. 

The  hyperbolic  fluxes  of  eqn(21)  are  discretized  using  a  Roe-type  approximate  Rie¬ 
mann  solver.  [20]  In  this  method  the  overall  solution  is  built  upon  the  solutions  to 
the  Riemann  problem  defined  by  the  discontinuous  jump  in  the  solution  between  each 
pair  of  cells.  The  numerical  flux  for  a  first-order  accurate  (in  space)  Roe-type  solver  is 
written  in  symmetric  form  as 


^1+1/2  =  X  (^i-hl  +  Fi)  —  “  ^  Ik  (Qi+1  “  Qi)  \^k\'^k 


(27) 


where  Vk  is  the  right  eigenvector,  Xk  is  the  absolute  value  of  the  eigenvalue, 
and  Ik  is  the  k^^  left  eigenvector.  The  values  at  the  cell  interface  (i  -1- 1/2)  are  obtained 
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Figure  3:  Schematic  of  the  arbitrary  shaped  three-dimensional  finite  volume  cell  to  be  used 
by  the  two-fluid  algorithm. 

either  by  a  simple  average  or  by  forming  a  more  accurate  Roe-average  of  the  neighboring 
cells.  The  flux  calculated  as  above  is  normal  to  the  cell  interface  which  is  the  desired 
orientation  for  applying  the  divergence  theorem  in  a  finite  volume  method.  A  typical 
three-dimensional,  finite- volume  grid  cell  is  shown  in  Figure  3. 

The  solution  of  Maxwell’s  equations  [eqn(4)]  is  required  at  each  time.  For  the  prelim¬ 
inary,  one- dimensional  algorithm  developed  in  this  stage  of  the  project,  a  magnetostatic 
model  for  the  electromagnetic  fields  is  assumed.  The  electric  field  is  determined  by  solv¬ 
ing  Poisson’s  equation  with  the  charge  density  and  the  introduction  of  the  electrostatic 
potential,  0.  Maxwell’s  equations  then  become 

V  •  E  =  -VV  =  eim  ~  ne)/eo.  (28) 

The  applied  magnetic  field  is  constant.  An  assumption  of  quasineutrality  {ui  =  rig)  is 
not  necessary. 

3  Project  Implementation  and  Results 

3.1  Algorithm  Details 

The  two-fluid  algorithm  has  been  successfully  implemented  in  one- dimension  with  an 
electrostatic  field  model.  Ampere’s  law  was  modified  by  including  an  electrostatic  po¬ 
tential  which  improves  the  quality  of  the  solution  for  larger  Xd  {^d  ^  simulations 
by  generating  electric  fields  that  simultaneously  satisfy  Ampere’s  law  and  Poisson’s 
equation.  The  stiffness  of  the  two-fluid  equations  was  relaxed  by  implicitly  treating 
the  non-homogeneous  source  terms.  The  homogeneous  fluxes  were  calculated  using  an 
approximate  Riemann  solver  for  each  fluid.  The  non-homogeneous  sources  terms  were 
then  solved  implicitly  to  complete  the  two-fluid  solutions. 

3.1.1  Maxwell’s  Equations  with  Electrostatic  Potential  Correction 

Maxwell’s  equations  couple  the  electron  fluid  with  the  ion  fluid  through  the  non- 
dimensional  field  equations.  Ampere’s  law  describes  the  time  evolution  of  the  electric 
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(29) 


field  in  the  presence  of  currents  and  magnetic  fields. 

60-^  =V  X  B//Xa-J 

However,  the  discretized  equations  do  not  guarantee  the  calcualted  electric  field  satisfies 
Poisson^s  equation  [eqn(28)]  which  describes  the  electrostatic  electric  field  generated 
by  a  collection  of  charges.  This  problem  is  similar  to  the  numerical  generation  of 
magnetic  fields  that  are  not  divergence-free.  Using  the  technique  of  Munz  [21],  MaxwelPs 
equations  were  re-written  to  include  an  electrostatic  potential,  <^,  which  propagated  the 
error  in  eqn(28)  out  of  the  domain  at  the  speed  P.  The  modified  equations  became 

eofio^  =  VxB  +  V<j>-j  (30) 

ot 

and 

^-r'*V-E  =  r=*e(ni-ne)/e„.  (31) 

at 


3.1.2  Homogeneous  Fluxes 

The  two  fluid  system  just  described  consists  of  two  hyperbolic  systems,  electron  fluid  and 
ion  fluid,  which  are  coupled  through  Maxwell’s  equations  source  terms.  It  is  important 
to  note  that  the  homogeneous  part  of  the  electron  system  depends  only  on  the  electron 
fluid  variables,  and  the  homogeneous  part  of  the  ion  system  depends  only  on  the  ion 
variables.  As  a  result,  separate  Riemann  solvers  are  used  to  solve  the  homogeneous  parts 
of  the  electron  fluid  and  the  ion  fluid.  Each  of  the  homogeneous  hyperbolic  systems  is 
written  in  the  form 

^+V-f’h  =  0  (32) 

where  Q  is  the  vector  of  piece-wise-constant,  conserved  variables  and  is  the  tensor 
of  hyperbolic  fluxes,  as  defined  in  eqns(22)  and  (23). 

The  approximate  Riemann  solver  solves  the  one- dimensional  Riemann  problems 
across  each  cell  interface.  Consequently,  it  is  only  necessary  to  include  the  flux  per¬ 
pendicular  to  a  face  in  our  discretization.  The  discretized  homogeneous  hyperbolic 
equations  are  solved  as  presented  in  Sec.  2.4.2. 

The  modified  Maxwell’s  equations  are  solved  with  an  approximate  Riemann  solver 
where  the  conserved  vector  is 

'b; 

Q=  E,  (33) 

4> 


and  the  flux  vector  is 


(34) 


For  each  system  a  matrix  A  is  calculated  such  that  AE  =  At^Q.  For  Maxwell’s 
equations,  A  is  simply  the  flux  Jacobian  matrix  |§.  For  the  fluid  systems,  A  is  the  Roe 
matrix  based  on  the  flux  Jacobian.  [22]  The  eigenvalues,  eigenvectors  and  characteristic 
waves  are  calculated  from  the  A  matrix. 
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3.1.3  Implicit  Treatment  of  the  Source  Terms 

The  stiffness  of  the  two-fluid  algorithm  was  generated  by  the  large  value  of  the  right 
hand  side  of  eqn(21).  Therefore,  the  source  terms  are  updated  implicitly  at  which  point 
the  three  separate  systems  (electron  fluid,  ion  fluid,  and  electromagnetic  fields)  are 
treated  as  one  large  system.  The  full  discretized  equation  with  implicit  source  terms  is 


-  Q" 
At 


^(FT+y2-Fr-i/2)+r^^ 


(35) 


where  is  the  source  term  evaluated  at  the  n  H-  1  time  step.  The  source  term  is 

expanded  in  a  Taylor  series  to  first  order  and  the  resulting  Newton  iteration  is  written 
as 


/_i_  _  _  A 

\At  SQV  ^  At  y 


{FP+i/2-Fr-y2)+i’’^ 


(36) 


where  AQ  =  and  k  is  the  iteration  variable. 

After  the  numerical  fluxes  are  calculated  explicitly  as  described  in  the  previous 
section,  the  source  terms  are  solved  implicitly.  Within  each  iteration  the  linear  system, 
eqn(36),  is  solved  using  a  symmetric  Gauss-Seidel  method.  The  flux  Jacobian  in  eqn(36) 
is  recalculated  and  Q  is  updated  after  each  iteration  until  the  2-norm  of  AQ  has  dropped 
by  several  orders  of  magnitude,  after  approximately  three  iterations.  It  is  important 
to  note  that  is  local,  and  t/;  does  not  depend  on  the  value  of  Q  in  adjacent  cells. 
This  means  that  the  computational  work  required  to  solve  eqn(36)  increases  linearly 
with  the  number  of  cells  in  our  domain,  even  when  the  algorithm  is  extended  to  three 
dimensions. 

In  many  cases  the  different  solvers  can  be  updated  using  different  time  steps.  For 
example,  for  MHD  like  problems  the  ion  equations  can  be  ignored  for  several  time  steps 
so  a  simplified  system  consisting  only  of  Maxwell’s  equations  and  the  electron  equations 
is  used.  After  several  time  steps  the  full  system  is  used  and  the  ion  fluid  is  updated. 
The  choice  of  time  step  depends  on  the  respective  stiffness  of  the  source  terms  for  each 
solver  as  well  as  the  hyperbolic  wave  speeds. 


3.2  Coplanar  Riemann  Problem 

The  preliminary  algorithm  was  tested  for  its  ability  to  capture  shock  and  its  ability  to 
capture  plasma  phenomena  such  as  Langmuir  oscillations. 

The  coplanar  Riemann  (or  shock  tube)  problem  is  used  to  test  the  fundamental  be¬ 
havior  of  the  approximate  Riemann  solver.  As  expected  from  an  approximate  Riemann 
solver,  wave  structures  are  well  captured.  The  behavior  is  best  seen  when  the  effect  of 
the  electromagnetic  fields  is  eliminated  by  artifically  setting  the  electron  and  ion  charge 
to  zero  and  reducing  the  effect  of  the  large  mass  difference  by  setting  the  electron  mass 
to  20%  of  the  ion  mass,  rrie  —  0.2mi.  The  electron  and  ion  fluids  then  decouple  and 
behave  as  independent  fluids.  Results  are  shown  in  Figs.  4  and  5.  The  coplanar  Rie¬ 
mann  problem  is  initiated  with  normalized  electron  and  ion  number  densities  of  4.0  on 
the  right  hand  side  of  the  domain  and  1.0  on  the  left  hand  side  of  the  domain.  The 
normalized  electron  and  ion  temperatures  are  initialized  at  a  uniform  value  of  10.  The 
domain  is  discretized  into  256  volumes.  Fig.  4  shows  the  evolution  of  the  electron  and 
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Figure  4:  Electron  and  ion  number  densities  for  a  shock  tube  simulation.  The  shock  wave, 
contact  discontinuity,  and  rarefaction  wave  can  be  seen  in  each  fluid. 


Figure  5:  Electron  and  ion  current  densities  for  a  shock  tube  simulation. 
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Figure  6:  Time  history  of  the  electron  velocity  measured  at  a  fixed  spatial  location  demon¬ 
strating  Langmuir  plasma  oscillations. 


ion  number  densities  after  150  time  steps.  The  characteristic  shock  wave,  contact  dis¬ 
continuity,  and  rarefaction  wave  are  clearly  evident  in  both  the  electron  and  ion  fluid. 
The  more  massive  ion  fluid  has  slower  characteristic  wave  speeds.  The  sound  speed 
varies  like  the  squareroot  of  the  mass.  Fig.  5  shows  the  evolution  of  the  electron  and 
ion  current  densities  in  the  longitudinal  direction  after  150  time  steps.  When  the  ef¬ 
fect  of  the  electromagnetic  fields  is  included,  the  faster  electron  fluid  is  slowed  by  the 
slower,  massive  ion  fluid.  The  number  density  jumps  are  altered  as  the  flow  is  limited 
by  ambipolar  motion. 


3.3  Langmuir  Oscillations 

The  ability  to  capture  gas  dynamic  waves  as  shown  in  the  coplanar  Riemann  problem 
is  expected  from  an  approximate  Riemann  solver.  The  approximate  Riemann  solver 
for  two  fluid  plasma  model  must  also  capture  plasma  wave  behavior.  This  ability  is 
demonstrated  by  simulating  Langmuir  plasma  oscillations.  Langmuir  plasma  oscilla¬ 
tions  are  the  plasmas  response  to  charge  concentrations  created  by  perturbations  in  the 
plasma.  In  one  dimension,  the  equation  of  motion  for  electrons  reduces  to  a  second 
order  differential  equation  for  the  particle  position. 


dPx 


rieC^ 

- X  = 

TTIqEq 


(37) 


where  ujpe  is  the  electron  plasma  frequency.  If  the  plasma  pressure  is  negligible  and  the 
fluid  velocity  is  sufficiently  small,  the  two  fluid  plasma  model  can  produce  Langmuir 
plasma  oscillations.  A  sinusoidal  perturbation  is  initialized  to  the  electron  fluid. 
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Figure  7:  Spatial  structure  of  the  electron  density  and  velocity  at  single  time  during  a 
Langmuir  plasma  oscillation  simulation. 


The  results  are  shown  in  Figs.  6  and  7.  The  theoretical  plasma  frequency  as  deter¬ 
mined  from  the  previous  equation  is  =  0.1466  for  the  normalized  plasma  parameters 
used  in  the  simulation.  The  numerically  generated  frequency  of  Upe  =  0.1462  agrees 
with  the  theoretical  value  with  an  error  of  0.27%. 

Though  the  results  presented  are  from  a  simplified  version  of  the  proposed  algorithm, 
the  success  and  accuracy  of  the  simulations  demonstrate  the  ultimate  potential  for  a 
new  algorithm  based  on  the  two-fluid  plasma  model. 

3.4  Debye  shielding 

In  this  problem  sheet  of  charge  is  introduced  near  the  center  of  the  domain  to  an 
otherwise  neutral  plasma.  The  electrons  and  ions  respond  to  shield  out  the  electric  field. 
The  modified  Maxwell’s  equations  are  used  in  this  simulation  because  charge  separation 
is  the  primary  effect  for  this  problem  which  has  A/j/L  =  0.1,  so  it  is  important  to 
maintain  the  correct  divergence  of  the  electric  field.  Furthermore,  when  the  divergence 
correction  is  not  used,  the  electric  field  at  the  center  of  the  domain  decays  and  the 
information  about  the  static  charge  distribution  is  lost.  This  is  a  dynamic  simulation 
which  is  stopped  when  the  numerical  dissipation  brings  the  solution  to  a  steady  state. 

For  the  simulation  results  presented  in  Figs.  8  and  9,  1000  grid  points  are  used. 
Without  the  sheet  of  charge  the  number  densities  would  be  uniform.  Without  the 
charged  fluids  the  electric  field  would  be  uniform.  The  analytical  solution  is  plotted  in 
Fig.  9  for  comparison. 

The  analytical  solution  assumes  constant  temperature  with  =  Te.  The  number 
densities  are  assumed  to  be  perturbed  such  that  Ui  =  nioCl—^^)  and  rig  =  neo(l+^)  where 
riio  =  Tieo-  Thermal  conductivity  is  neglected  from  the  calculated  two-fluid  solution,  so 
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Figure  8:  Spatial  structure  of  the  electron  and  ion  number  densities  showing  the  charged 
fluids  response  to  the  sheet  of  charge. 


Figure  9:  Spatial  structure  of  the  electric  field  showing  the  shielding  effect.  Without  the 
charged  fluids  the  electric  field  would  be  uniform.  The  analytical  solution  for  this  problem 
is  shown  by  the  dotted  fine  for  comparison. 
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Figure  10:  The  variation  of  the  Child-Langmuir  space  charged  limited  current  as  the  dis¬ 
tance  between  the  electrodes  are  adjusted.  The  data  points  are  from  numerical  simulations, 
and  the  solid  curve  is  predicted  by  the  theory  [eqn(38)]. 


the  steady  state  solution  is  not  in  perfect  thermal  equilibrium  which  may  explain  the 
discrepancy  with  the  exact  solution  in  Fig.  9. 


3.5  Child-Langmuir  Current  Saturation 

The  most  significant  accomplishment  during  this  grant  period  has  been  the  development 
of  the  approximate  Riemann  solver  for  the  two-fluid  plasma  model.  This  section  briefly 
describes  the  progression  of  the  development. 

A  fundamental  limitation  of  MHD  algorithms  is  the  inability  to  capture  inherently 
two  fluid  behavior,  such  as  space  charge  effects.  Space  charge  effects  can  occur  when 
finite  electric  potentials  are  applied  to  the  plasma  or  when  the  plasma  comes  in  contact 
with  insulating  or  conducting  walls.  Space  charge  effects  can  be  investigated  by  consid¬ 
ering  an  applied  electric  potential  across  a  vacuum  with  a  plasma  fluid  at  the  boundary. 
The  electron  plasma  fluid  (for  example)  is  extracted  from  the  boundary  and  accelerated 
through  the  potential  until  it  leaves  the  opposite  boundary  of  the  domain.  Since  the 
plasma  is  extracted  from  the  boundary  with  negligible  velocity  before  it  is  accelerated, 
the  electron  fluid  density  increases  and  a  space  charge,  virtual  cathode  tends  to  form 
at  the  boundary.  The  virtual  cathode  retards  the  extraction  of  additional  plasma  from 
the  boundary.  The  space  charge  effect  limits  the  extracted  current  according  to  the 
Child-Langmuir  law. 


3  CL 


V  ^ 


(38) 
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which  determines  the  saturation  current  obtained  for  an  electrode  spacing  of  L  and  an 
applied  potential  of  (j). 

Child-Langmuir  current  saturation  is  simulated  with  the  two  fluid  algorithm.  The 
results  are  presented  in  Fig.  10.  The  data  presented  illustrate  the  variation  of  the 
saturation  current  as  the  spacing  between  the  electrodes  is  adjusted.  The  data  points 
are  generated  from  the  numerical  simulations,  and  the  solid  curve  is  the  plot  of  eqn(38). 
The  error  in  the  fit  at  the  extreme  values  of  the  potential  are  the  result  of  numerical 
dissipation  of  the  first-order  flux  calculation.  It  is  expected  that  the  fit  would  improve 
with  a  higher-order  method. 


4  Professional  Interactions 
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Thomas  R.  Jarboe  Professor 

Chris  Aberle  Graduate  Student 

John  Loverich  Graduate  Student 

4.2  Collaborations 

4.2.1  Air  Force  Research  Laboratory 

Dr.  Robert  Peterkin,  Jr.  of  the  Electromagnetic  Sources  Division  of  the  Air  Force 
Research  Laboratory  at  Kirtland  AFB  on  three-dimensional  multigrid  algorithms  for 
MACH3  an  development  of  a  parallel  PIC  (particle  in  cell)  code  for  microwave  simu¬ 
lations.  Knowledge  developed  in  the  area  of  relaxation  schemes  was  implemented  into 
the  ICEPIC  code  to  make  a  3-D  Poisson  solver.  The  solver  was  needed  to  determine 
electric  field  concentration  on  a  high  power  microwave  source. 

4.2.2  University  of  Washington 

Prof.  Scott  Eberhardt  of  the  Aeronautics  and  Astronautics  Department  and  Prof.  Randy 
LeVeque  of  the  Applied  Math  Department  on  approximate  Riemann  solvers  and  their 
applications  to  multidimensional  problems,  in  particular,  the  correction  to  transverse 
fluxes. 

Prof.  Tom  Jarboe  of  the  Aeronautics  and  Astronautics  Department  on  the  MHD 
stability  simulations  of  spheromaks  in  realistic  three-dimensional  geometries  and  the 
beta  pressure  limit  of  driven  spheromak  plasmas.  This  is  an  ongoing  collaboration  that 
resulted  in  several  publications  and  an  experimental  project  for  Prof,  Jarboe. 


4.3  Publications 

Two  papers  describing  this  project  and  related  research  were  presented  at  the  AIAA 
Computational  Fluid  Dynamics  conference.  One  paper  titled  “An  Approximate  Rie¬ 
mann  Solver  for  MHD  Computations  on  Parallel  Architectures”  provided  an  overview 
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of  the  project.  The  other  paper  titled  “Application  of  Analytical  Methods  to  Com¬ 
puting  Numerical  Flux  Jacobians”  described  the  work  on  the  analytical  flux  Jacobian 
calculation. 


5  Conclusions 

The  research  accomplishments  of  this  research  effort  represent  significant  advances  in 
the  field  of  plasma  dynamic  simulations.  This  advancement  is  important  because  plasma 
researchers  are  increasingly  finding  that  the  assumptions  made  in  the  single-fluid  MHD 
model  are  too  restrictive  to  model  selected  plasma  dynamic  phenomena,  and  currently 
two-fluid  plasma  simulation  codes  have  been  limited  to  linearized  models  with  only  two 
spatial  dimensions.  The  accomplishments  indicate  that  the  proposed  research  goals  are 
being  met. 

The  successful  development  of  a  one-dimensional  algorithm  based  on  the  two-fluid 
plasma  model  demonstrates  the  potential  of  the  project,  and  it  indicates  that  this 
project  is  reaching  its  objectives.  The  research  related  to  this  project  has  been  presented 
at  international  conferences.  Valuable  collaborations  have  been  formed  with  the  Air 
Force  Research  Laboratory. 

The  continuing  development  of  this  project  will  include  making  the  algorithm  more 
robust  to  handle  electromagnetic  fields,  develop  an  electromagnetic  solver  for  arbitrary, 
cell-centered  grids,  and  extending  the  algorithm  to  three-dimensions. 
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